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TECHNICAL PUBLICATION 


SYSTEM SIMULATION BY RECURSIVE FEEDBACK: COUPLING A SET 
OF STAND-ALONE SUBSYSTEM SIMULATIONS 

1. INTRODUCTION 


In design and development of space systems, digital dynamic simulations have potential to reduce 
physical testing requirements, increase confidence in systems that cannot be tested, and ensure that sub- 
system requirements are not more stringent than necessary. However, simulations that address interactions 
of multiple disciplines and subsystems can be difficult to build. Knowledge and expertise required at the 
physics level, especially for high-fidelity developments, do not normally reside within one person or orga- 
nization for more than one discipline or subsystem. Simulations involving individual disciplines or sub- 
systems usually exist or can be readily constructed, and system engineers often understand the qualitative 
nature of interactions among them; but a means of obtaining quantitative-coupled system-level results is 
often unavailable. Therefore, specific practical approaches to construction of high-fidelity system-level 
simulations are being developed from numerous perspectives across the engineering community. 

Digital dynamic simulations are often built by deriving or constructing differential -algebraic equa- 
tions that define the system model and implementing algorithms to solve those equations numerically, as a 
single coupled set, using one of several standard or modified numerical integration methods. In some 
cases, software is written in a language such as FORTRAN or C to define the equations, and an existing or 
slightly modified ordinary differential equation (ODE) solver, perhaps a Runge-Kutta implementation, is 
employed to perform the integration. In other cases, an in-house development such as Marshall System for 
Aerospace Simulation (MARS YAS ) or a commercial off-the-shelf product such as MATLAB®/SIMULINK® 
is used to act as a higher level facilitator of what is ultimately a mathematically similar approach. 

System-level digital solution of coupled “stand-alone” dynamical simulations is a fundamentally 
different approach to system simulation, and fundamentally different numerical procedures are required. 
Recursive feedback is a conceptual method from which a family of appropriate algorithms, processes, and 
specific numerical procedures can be derived. 

Coupled-subsystem-simulation architecture provides near complete independence of stand-alone 
subsimulations and naturally facilitates high fidelity and broad scope through collaboration across inter- 
faces that can be implemented in the same physical and engineering terms that define them in the real 
system. Such simplicity promotes clarity of communication and ease of understanding, both of which have 
many positive benefits. Individual subsimulations can potentially be implemented on separate, remote, and 
dissimilar computational platforms, and this portends numerous advantages and possibilities as computer 
network capabilities improve. 



2. RECURSIVE FEEDBACK ALGORITHM 


2.1 Background 

Brief discussions of high-fidelity, multidiscipline simulation are set forth in references I and 2. 
which resulted fiom association with the Structures and Dynamics Laboratory High-Fidelity Multidisciphne 
Simulation Development team. Reference 1 attempts to illustrate the long-term goal of multidiscipline 
simulation by example, while reference 2 presents the goal in a more generalized way. Reference 3 docu- 
ments development ot an earlier FORTRAN general purpose dynamical simulation tool (ANASIM) that 
uses recursive feedback as the system solution algorithm. This tool allows one to construct a simulation by 
choosing system elements from a menu and specifying connections (signal flow) among them. The ele- 
ments (integrators, summers, limiters, etc.) are, in effect, small stand-alone subsimulations. A definition 
and explanation of the method is presented, and an analytical demonstration of response computation for a 
first-order, single-state, constant-coefficient system with a step input is included. The demonstration, or 
mathematical “show that,” is accomplished through symbolic algebraic application of the algorithm fol- 
lowed by direct comparison of results with the exact closed-form solution. 


In reference 4, the method is analytically demonstrated for some simple examples within a variety 
of system classes. The first is the /(-dimensional linear constant-coefficient case defined by the matrix 
vector equation .\{t) = A.\(t ) + Bu(t ) , which is a standard state-variable form. The remaining four are sca- 
lar (single state) cases that include two nonlinear, continuous cases; a linear time-varying coefficient case; 
and a numerically demonstrated case involving a state-dependent discontinuity and an algebraic loop. The 
analyses of reference 4 provide additional mathematical insights into the process, particularly with respect 
to convergence. 

Analysis of the method is extended to coupling of stand-alone subsystems in references 5 and 6. 
Reference 5 mathematically defines and analytically demonstrates recursive feedback for generation of 
responses from a generic linear-dynamic system composed of internally inaccessible subsystems, and ref- 
erence 6 continues beyond to include subsystems that have potential to produce algebraic loops at the 
system level because they contain direct feed-through of input. 

References 5 and 6 specifically address process convergence and essentially guarantee that recur- 
sive feedback can be used as a system-level solution algorithm for simulation of linear systems composed 
of coupled stand-alone subsystems for a broad range of ordinary circumstances. Combined with the analy- 
ses of reference 4 and experience associated with ANASIM (ref. 3), it seems reasonable to expect that the 
method can be dependably applied to nonlinear system simulation. Only experience, and perhaps more 
analysis, will determine the effectiveness of extrapolation to these supposedly more stringent conditions. 





2.2 Concept 


Referring to figure I . the subsystems are represented by independent, self-contained, time-domain, 
dynamical simulations that can be commanded to run from a given set of initial conditions over a predeter- 
mined segment of time when provided with input signals as functions of time over that interval. In general, 
the time segment, also referred to as a convergence interval, is short compared with the total desired system 
simulation duration, but may be significantly larger than the integration step size associated with a typical 
subsimulation. For each new stage of recursion, subsystem input signals are computed using original 
system-level input signals summed with system-level feedback (coupling) signals generated at the 
current stage. New (next recursion stage) subsystem responses are computed using new subsystem 
input signals in conjunction with the same (original) initial conditions. The process continues recursively 
until convergence is achieved. After convergence, initial conditions are replaced with final conditions and 
the process is restarted for the next segment of time. System response over the total desired duration of the 
simulation is the concatenated set of combined subsystem responses over many time segments. 


Subsystem Inputs Subsystem Responses 



Figure 1 . Schematic of recursive feedback process. 
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2.3 Insights for Specific Design 


Because figure I describes a system of considerable generality whose complexity for a specific 
case can range from simple to great, further clarification of the process seems appropriate. The simplest 
“system” could arguably be a single integrator in the forward path, a constant multiplier (gain) in the 
backward (feedback) path, no (zero) system-level input, and a nonzero initial condition on the output. In 
that case, there is one linear subsystem, and that subsystem is a single integrator. The output function is the 
integral of the subsystem-level input function with respect to time, and can be generated independently and 
directly by simply computing the area under the input curve over the time interval of interest (convergence 
interval) and adding the initial condition. The coupling matrix becomes a scalar (value of the loop gain), 
and actually lepresents direct feedback rather than coupling. The process, however, does not explicitly 
make such a distinction, nor does it need to. The summing junction becomes moot since there is no system- 
level input. The recursive transfer point remains as depicted in figure 1, but does not exist in the analog 
system to be digitally simulated. The scalar differential equation that models the analog system is 
)’(/) — k}’( t ) , where t is time and k is the loop gain, and the exact closed-form solution for the output 
response is given by the exponential function 


y(D = y(0) \e k '] 


( 1 ) 


where v(0) represents the initial value of the output v(/). 


Applying recursive feedback to this system, the output is first computed as a function of time over 
the convergence interval as if there were no (zero) feedback. This function is referred to as the stage zero 
output response. Because the system input function is also zero in this case, the response is a constant 
function of time symbolically represented by 


y' o, (0 = y(0) , (2) 

where the parenthesized superscript is the stage index. Next, the stage zero feedback function is computed 
by multiplying the stage zero output response function by the loop gain k. After “moving” the feedback 
function data across the recursive transfer point and summing it with the system-level input function (zero), 
it becomes the stage one subsystem input function. Stage one response can now be computed by numerical 
integration (area under the curve) of the new (stage one) input function while applying the original initial 
condition v(0). A symbolic representation of the result is 

y" > (0 = >’(0)[l + £r] . (3) 

Repeating the process, symbolic representation of the stage two result is 


y 2) (r)= v(0)[l + kt + 


an 2 

2! J 


(4) 


4 



By induction, the nlh stage response is given by 


y'"'(t) = y(0) 


\ + kt + 


(MY 


( kt )" 

+ ■■■ + — 

nl 


(5) 


Clearly, if sufficiently continued, the recursive process yields a solution that approaches the exact 
solution >’(/) = y(0) [e kl ] . In theory, the time interval (convergence interval) for this case can be as long as 
desired; in practice, it must be restricted because of numerical integration and roundoff errors. As men- 
tioned earlier, system response over the problem duration is represented by a concatenation of responses 
over as many convergence intervals as needed, and a single convergence interval is measured in multiple 
integration steps. It is noted that each recursion stage improves accuracy of the approximate system 
response by adding the next higher order term in a series expansion of the exact system response. This type 
of behavior is seen in all of the linear system analyses of references 3-6. For the nonlinear systems of 
reference 4, each recursion stage adds higher order terms; but for a given number of correct terms in the 
series expansion, a greater number of recursion stages is required. For these systems, the convergence 
interval cannot be arbitrarily long, even without integration and roundoff errors, and the degree of restric- 
tion required depends on the initial condition. 

Considering the system in figure I in its more general context, the off-diagonal elements of the 
coupling matrix represent the means by which output from one subsystem is processed and channeled to 
construct input for another. Because direct dependence of the input to any subsystem on the output of that 
same subsystem is likely to have been handled internally, the elements on the diagonal are likely to be zero; 
although, as demonstrated, the method does not require this to be the case. In development of a system 
simulation by coupling of subsystem simulations, the nature of the coupling matrix is, to a large degree, a 
design choice. 

Consider, for example, the “system" to be an undamped oscillator governed by the scalar differen- 
tial equation y(t) + co 2 y(t) = 0 , where a> is the frequency of oscillation. A simulation of this system by 
recursive feedback can be set up in at least two very different ways. One design possibility is to define a 
single subsystem, best described as a double integrator, or two single integrators in series. The output 
function of the subsystem is determined by integrating its input function twice. The coupling matrix once 
again degenerates to a scalar constant, -<y 2 , and the system input function is zero. Assuming a nonzero 
initial condition y(0) on the output (displacement), application of recursive feedback to this system results 
in a truncated series definition of y(t), and at the «th recursion stage, a symbolic representation is 


y'"’(')=.v(0) '-^ + ^-’- + (-i)' w ''^ 


( 6 ) 


Then for a sufficiently large number of recursion stages, the approximate solution approaches the 
exact solution y(t) — v(0) [coscur] . Once again the number of terms in the series corresponds directly to 
the number of recursion stages, but the order of the series with respect to time is twice the number of 
recursion stages. 
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In contrast, a second design possibility is to define two subsystems that are both single integrators. 
The governing differential equation is now given by the two-dimensional matrix-vector equation 


Y{t) = KY{t) , 


( 7 ) 


the coupling matrix is given by 


and the output function is given by 




(8) 


Y(l) = 


v,(0 

y 2 (t)\ ' 


( 9 ) 


where y j (t) represents displacement and represents velocity. Application of recursive feedback to this 

system results in a truncated series definition of F(r). At the mh recursion stage, the numerically computed 
approximate system response function is symbolically represented by 


Y { "'(t) 


l+Kt + K : — + 

21 


+ K"~ 


Y ( 0 ) . 


( 10 ) 


Again, each numerical recursion stage, in effect, adds one term to the truncated series approxima- 
tion of the exact solution. 


For the two-subsystem, twin-integrator approach, it takes two recursion stages to reach the same 
order of approximation in time that was reached in one stage with the single-subsystem, double-integrator 
approach. However, the twin integrators represent independent operations that can be performed in paral- 
lel, while the double integrator represents two operations that must be performed in series. In the double 
integrator case, velocity, y(t), does not appear at the system level, only one function crosses the recursive 
transfer point, and only one function is checked for convergence. In the twin integrator case, velocity, y,(t), 
appears at the system level, both Vj(/) and y 2 (0 cross the recursive transfer point, and both functions can be 
checked for convergence at the system level. Clearly, in a more complex situation, the relative merits of 
these approaches become simulation design and performance issues. 

For the broader range of systems contained within the framework of figure 1, the nature of the 
coupling matrix is a function of many possible simulation design choices. In its most basic form, a coupler 
would simply route signals from one subsystem to another without modifying them. In its most complex 
form, it could become a multidimensioned nonlinear operator as well. Typical coupler functions are likely 
to be coordinate transformations, gain multiplications, interpolations, and other data modification or rout- 
ing tasks that must be performed to make proper interface connections among a set of predefined or exist- 
ing subsimulations. 
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As implied by the double integrator example (when viewed as two single integrators in series), it is 
also possible to distribute the coupling functions throughout the system. In a particular situation, one might 
choose to view several subsimulations as a combined set if the nature of the coupling among them is such 
that an output response can be determined from an input function without necessity for recursion (there is 
no feedback loop internal to the set). In that case, processing information as it is transferred from one set 
member to another represents a coupling function that does not reside as an element of a centralized cou- 
pling matrix. Utilizing the powers of block diagram algebra or signal flow analysis, many arrangements ot 
subsystems and couplers may be possible for representation of a given system. Such design freedom can 
have a significant impact on the configuration, performance, and efficiency of the simulation. 

To be of practical value in dynamical systems simulation, the recursive feedback concept must be 
applicable to nonlinear models. The “guarantees’' and insights that come with analyses like that of refer- 
ence 6 are weakened at best and totally lost at worst without the linear model assumption. However, one 
does have a guarantee of solution existence and uniqueness for a broad range of sufficiently well-behaved 
nonlinear systems (ref. 7). A converged recursive feedback process is indicative ot a solution with excep- 
tion of questions concerning discretization errors, roundoff errors, inappropriate tolerances, etc. The 
potential for convergence can normally be improved by reducing the length of the convergence interval, 
not unlike reducing integration stepsize in a conventional situation when satisfactory results are not achieved. 
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3. EXAMPLE SYSTEM 


The example system is designed to be simple but to reflect real possibilities. It consists of a satellite 
that has an attitude control system, is in low-Earth orbit, and is experiencing aerodynamic effects in addi- 
tion to gravitational effects. The spacecraft is idealized as a rigid body of rectangular “box” shape, illus- 
trated in figure 2. The center of mass is offset from the geometric center by a small amount in all three axes, 
and moments of inertia are such that there are two large but somewhat unequal principal values, while the 
third principal value is much smaller. The principal axes are normal to the box surfaces so that all cross 
products of inertia are zero. The models are based on constant density atmosphere, spherical Earth (inverse 
square) gravity, and linear attitude control for small motion about an inertial reference. Attitude control can 
be partially or totally removed so that the body is allowed to tumble under the influence of gravity and 
aerodynamics. Gravitational, aerodynamic, and control torques affect the rotational dynamics of the rigid 
body, while gravitational and aerodynamic forces affect the translational dynamics. Both translational and 
rotational dynamics affect aerodynamic forces and torques, so a feedback loop is apparent. This example 
falls in the general category of nonlinear, continuous systems. 


Aerodynamic Data: 

Drag Coefficient: 2.25 
Density: 1.95 xIO -11 kg/m 3 


Dimensions: 

Length ( L ): 15 m 
Width (W): 4 m 
Height (W): 3 m 


Center of Mass Offset: 

X=0.1 m 
Y = 0.1 m 
Z- 0.1 m 


Notes: 

1. © depicts center of mass location. 

2. Roll, pitch, and yaw are Euler angles about the 
X, V, and Z axes, respectively. The sequence 
is yaw -pitch -roll. 

3. All axes penetrate facets at their geometric centers. 


Inertial Pointing Control Gains: 

Position: 1,000.0 N-m/rad, All Axes 
Rate: 100.0 N-m-sec/rad. All Axes 



Mass Properties: 

Mass: 10,000 kg 
Moments of Inertia: 
11,250.0 0 

0 193,125.0 

0 0 


0 

0 

173,812.5 


kg - m 2 


Figure 2. Spacecraft properties and geometric body-fixed coordinate system. 
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4. MATHEMATICAL MODELS 


4.1 Orbit Dynamics 


The orbit dynamics is modeled by 


where 


d 2 R 


ni 


~ F(»rav "b Faeru ' 


sat , 1 * grav ^ 1 aero 

dt 


p _ Earth w sat p 

r grav - _3 K ‘ 


(ID 


(12) 


G is the universal gravitational constant, and »? sat are masses of the Earth and the satellite. R, the 

orbit radius vector, extends from the center of the Earth to the center of mass of the satellite and is 
expressed in an Earth-centered reference frame that is defined to be nonrotating and assumed to be 
nonaccelerating. F aero , the aerodynamic force, is derived from an aerodynamics model to be discussed 
later. This development, as well as the attitude dynamics development, is standard and of the general type 
found in many dynamics texts such as reference 8. 


4.2 Attitude Dynamics 


The rotational dynamics is modeled as a rigid body by 


1(0 + (o x lo) — T contro] + Tg„ + T aero , 

where I is the moment of inertia matrix, d) is the angular velocity vector, T gg , is gravity gradient torque, 
Taero * s aerodynamic torque, and all quantities are expressed in a body-fixed coordinate system whose 
origin is at the center of mass. 


The attitude is defined by a set of Euler angles [$ v ,0 v ,</>- that indicate angles traversed by a hypo- 
thetical reference frame as it is rotated from coincidence with an inertial (nonrotating) attitude reference 
frame to coincidence with the body frame using three single-axis rotations in a r-y-.v sequence. The instan- 
taneous Euler angle rates are then related to body rates by 


</>v 


1 

sin (j) v tan <p v cos <p y tan (p v 

(Ox 

4 

= 

0 

cos <p x -sin (f) x 

(Oy 

4 


0 

sin0 v sec0 v cos</> v sec 0 V 

( 0 - 
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Although this particular formulation has a singularity at (f> x — iTT /2, it is valid except at that point, 
and can serve as the basis for a model if an adequate neighborhood of that point is avoided. 

A vector expressed in inertial coordinates is transformed to body coordinates through multiplica- 
tion of the corresponding column of inertial components by the direction cosine matrix 



1 

0 

0 

cost/> v 

0 

-sin0 v 

CO S<p- 

sin0. 0 

A = 

0 

cos 0 V 

sin0 v 

0 

1 

0 

-sin0- 

COS 0- 0 


0 

-sin0 v 

cos0 v 

sin0 v 

0 

COS0 V 

0 

0 1 


The attitude control system is modeled by 


1 control 


fl oA +fl iA 

C/()y0V + U\ yCOy , 

a Q: (l) : +a\ : (D : 


(16) 


where <p, to, and a are associated with attitude angles, body rates, and control gains, respectively, and the 
subscripts are self-explanatory. 


The gravity gradient torque is modeled by 9 


T = 

*2{Z 


Gm 


Earth 


R 


r x / ■ r 


(17) 


where r is a unit vector corresponding to R expressed in body coordinates. 

4.3 Aerodynamics 

A box with six rectangular side panels (facets) defines the satellite shape (fig. 2). Aerodynamic 
drag forces associated with the /th panel are approximated empirically by the formula 9 


Fi=< 


--pV^CrfA/inj -v)v 


0 


if h f > 0 

if h ■ v < 0 


(18) 
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where p is atmospheric density, V is velocity magnitude relative to the atmosphere, Cj is a drag coefficient, 
A is panel area, n is a unit vector normal to the panel and directed outward from the box, and v is a unit 
vector corresponding to V. If rotation of the Earth is neglected, no distinction between relative and abso- 
lute velocity need be made. The total aerodynamic force is given by 

F=ltl f 5' (l9) 

While the vector relationships of equations ( 1 8) and ( 19) are independent of any particular coordi- 
nate system, the most convenient choice is a body-fixed “geometric” coordinate system (fig. 2) whose 
origin and axes are defined from geometry (shape) considerations alone. The aerodynamic force as used in 
equation (1 1) is most appropriately expressed in inertial coordinates. Appropriate coordinate transforma- 
tion depends on the matrix A. defined by equation (15), and on the orientation of the body frame with 
respect to the geometric frame. As illustrated in figure 3, the axes of the body frame are not rotated with 
respect to the geometric frame, so the relevant transformation is identity. 



Figure 3. Relationship of body frame to geometric frame. 


1 I 



The total aerodynamic torque about the origin of the geometric frame is calculated from six con- 
tributing components and is given by 

T = Xti r cp i x F, . (20) 

where r cpi is a vector from the origin of the geometric coordinate system to the center of pressure for the /th 

facet. The total aerodynamic torque about the center of mass of the spacecraf t, as used in equation ( 1 3), is 
given by 

Tcn, = T — (f cm xF) . (21) 

where r cm is a vector from the origin of the geometric coordinate system to the center of mass (fig. 3). 
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5. SUBSYSTEM SIMULATIONS 


Simulations for each of the three subsystems (orbit dynamics, attitude dynamics, and aerodynam- 
ics) have been constructed as functionally separate entities, although they are implemented as parts of one 
single FORTRAN computer program. The attitude and orbit dynamics simulations both have 3 degrees of 
freedom, and each has its own separate implementation of a Runge-Kutta numerical integration process. In 
both cases, the three second-order differential equations have been rearranged to six first-order equations 
for interface with the integration algorithm. The aerodynamics simulation has no integration process 
because the model is purely logical-algebraic. All three simulations are “stand-alone" because each, in 
principle, is capable of producing output signals (functions of time) from input signals and initial condi- 
tions. Initial conditions, of course, do not apply to the logical-algebraic aerodynamic simulation. Input and 
output signals are represented by a sequence of linearly connected points equally spaced in time. Addi- 
tional points are obtained through interpolation by each simulation if needed. 

The orbit dynamics simulation input signal is aerodynamic force, and output signals are orbit 
radius and velocity vectors. The model, defined by equations (11) and (12), considers a point mass, inverse 
square gravity, and an external force (aerodynamics) that is an arbitrary but known function of time. All 
orbit dynamics simulation variables are referenced to an Earth-centered inertial coordinate system. 

The attitude dynamics model input signals are the orbit radius vector expressed in inertial coordi- 
nates, and aerodynamic torque expressed in body coordinates. Output signals are body rates and attitude 
angles. The model, defined by equations (13 ) — ( 1 7), considers rotational motion of a rigid body, gravity 
gradient torque, attitude control system torque, and an external torque (aerodynamics) that is an arbitrary 
but known function of time. 

The aerodynamics simulation input signals are orbit radius vector in inertial coordinates and orbit 
velocity vector in body coordinates. Output signals are aerodynamic force and torque expressed in geomet- 
ric coordinates. The model is logical-algebraic and is defined by equations (18)— (21 ). 
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6. SYSTEM SIMULATIONS 


Three system simulations have been constructed. The first simulation uses a centralized coupling 
implementation of recursive feedback, while the second simulation uses a distributed coupling implemen- 
tation. The third simulation uses a conventional first-order, single-equation set approach and is the standard 
against which the recursive feedback approaches are compared. 

By design, the function of a coupler is to perform coordinate transformations or other types of data 
conversions needed to make consistent, correct connections among the subsystem simulations. A single 
coupler can peiform all coupling functions from one central location, or a set of couplers that perform less 
extensive functions can be distributed throughout the system. Coupling signals, like input and output sig- 
nals, are represented by a sequence of points equally spaced in time. 


Figures 4 and 5 illustrate the centrally coupled system simulation. The subsimulations represent 
independent processes that can, in principle, be executed in parallel over the timespan of each convergence 
interval. Because the primary purpose is to investigate functionality rather than efficiency, no attempt has 
been made to use more than one processor in the single-computer FORTRAN implementation. Therefore, 
the parallel computation advantage is lost, and quantitative comparison of efficiency with the other simula- 
tions is not possible. 

Figure 6 illustrates the distributively coupled system simulation. Before applying the idea of recur- 
sive feedback, the subsystems have been arranged so that, to the extent possible, feedback is minimized. 
Because there is an information loop that ultimately must be closed, it is not possible to completely elimi- 
nate feedback. The subsimulations are identical to those of the centrally coupled system, but the system- 
level numerical process is substantially affected because, among other things, the number of scalar signal 
paths that cross the recursive transfer point has changed. This simulation was developed in part because the 
centrally coupled simulation, developed first, is noticeably slower than the standard method benchmark. 
More importantly, it serves to confirm functionality of the different coupling architectures. 

The conventional simulation is implemented from the combined set of equations that define the 
subsystem models. They have been collected and rearranged, by hand, to a single set of 12 first-order 
differential equations, an arrangement sometimes called a “state variables” approach to system formula- 
tion. In that context, the state vector is defined as 


X = [R\V \(l>\ 0)] T , (22) 

where /?, V 7 , <p , and u> are three-element row vectors representing orbit radius, velocity, Euler angles, and 
body rates. The system differential equations can now be represented as a first-order set by 


X = f(X, t) , 


(23) 
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Figure 4. Centrally coupled simulation diagram. 
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*** Geometric Coordinates 


Figure 6, Distribut i vely coupled simulation diagram. 

where time has been explicitly included only for completeness of the form. In the current simulation, time 
is only implicitly involved, because there are no system forces or torques that are explicit functions of time. 
A fourth-order Runge-Kutta numerical integration algorithm specifically designed to solve systems in the 
form of equation (13) is employed to propagate the system responses. 

All system and subsystem simulations were specially developed for this project. They are written in 
FORTRAN, all data were generated on the TD55/eds2 UNIX platform at Marshall Space Flight Center 
(MSFC). Plots were originally produced with an in-house program (developed by Warren Adams, TD55, 
MSFC) on the same machine. Plot images were migrated to a Windows desktop PC for publication. A full 
complement of plots is included in appendix A, and listings of the simulation programs are included in 
appendix B. 
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7. NUMERICAL RESULTS 


7.1 Run Case Definitions 


All data cases simulate motion of the same spacecraft starting from the same orbit and attitude 
related initial conditions. The orbit is initialized so that in the absence of aerodynamic drag it would be 
circular at an altitude of 300 km and an inclination of 30 deg. The initialization point is in the X-Z plane ot 
the coordinate system shown in figure 7. Spacecraft mass properties, dimensions, and control gains are 
given in figure 2. Initial attitude angles and body rates are zero. Other data constants and physical constants 
are identified in the computer listings. Case variations are developed from three attitude control system 
configurations and three simulation methods resulting in a nine-element run matrix. The attitude control 
configurations are listed as follows: 

(1) Full inertial control 

(2) Partial inertial (two-axis) control 

(3) No control (tumble). 

Case numbers correspond directly to the control system configuration list entry numbers, and method 
variations are associated with each case. The simulation methods are: 

(1) Centrally coupled recursive feedback (RF-C) 

(2) Distributively coupled recursive feedback (RF-D) 

(3) Conventional (unified) Runge-Kutta (RK-U). 


t Y tz 



Figure 7. Inertial coordinate system/initial orbit plane geometry. 
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Numerical algorithm parameter settings were the same for all cases. Both recursive methods used 
five sample points per convergence interval (including both end points), and the interval length was 0.4 sec. 
Relative (normalized) and absolute convergence tolerances were both 10~ 12 for the RF-C and 10~ 14 for the 
RF-D. The orbit and attitude subsimulations both used a fourth-order Runge-Kutta algorithm with a step 
size ot 0.1 sec, while the aerodynamics subsimulation computed aerodynamic effects at a 0. 1 -sec sample 
interval. The conventional simulation used the same Runge-Kutta algorithm with a step size of 0.1 sec. 

Results tor all possible (nine) combinations of control system configuration and solution methods 
were generated. For case I, in the absence of any well-defined requirement, a set of control system gains 
was chosen that produce reasonable (small angle) attitude dynamics responses. For case 2, the Z-axis 
control gains were set to zero; and for case 3, all control gains were set to zero. Simulation responses have 
been generated for 60,000 sec of problem time using all solution approaches for all cases. 

7.2 Response Comparisons 

The responses for case I (full control ) are identical over the full 60,000-sec duration of the simula- 
tion runs for all solution approaches (RF-C, RF-D, and RK-U). Appendix A.l contains 108 plots that 
compare RF-D results to conventional RK— U results for case 1. Radius, velocity, attitude angles, body 
rates, aerodynamic forces, aerodynamic torques, and altitude are shown as functions of time. The general 
format ot the appendix presentation is to show each response over time intervals of 0-20,000, 20,000- 
40.000, and 40,000-60,000 sec, and additionally, to show some selected plots that cover smaller time 
intervals of particular interest. 


The responses for case 2 (two-axis control) are, like case 1 , identical over the full 60,000-sec dura- 
tion of the simulation runs for all solution approaches. Appendix A.2 contains 92 plots that compare RF-D 
results to conventional RK-U results for case 2. As in the appendix presentation for case 1 , radius, velocity, 
attitude angles, body rates, aerodynamic forces, aerodynamic torques, and altitude are shown as functions 
of time, and the presentation format is also similar to that used for case I . 

For case 3 (tumble), the RF-C, RF-D, and RK-U simulation responses all agree for approximately 
45,000 sec of the duration; after that, the RF— C and RF— D responses begin to diverge from the RK— U while 
they continue to agree with each other for the full 60,000 sec. Figures 8-16 present plots selected specifi- 
cally to illustrate the early agreement and subsequent deviation of responses generated recursively from 
those generated conventionally. Attitude angles are shown in figures 8-10, body rates are shown in figures 
1 1-13, and altitude is shown in figures 14-16. Only the 40,000- to 60,000-sec time period is covered 
because no earlier deviation is apparent. Appendix A. 3 contains 60 plots that compare results of the RF— D 
approach to the conventional RK-U approach for this case. 
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Figure 8. Case 3. Euler angle, X axis (deg). 
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Figure 9. Case 3, Euler angle, Y axis (deg). 
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Figure 16. Case 3, altitude, third time segment (km). 


Altitude results for all control system configurations and two solution methods (RF-D and RK-U) 
are superposed in figures 17-19. The main point of illustration is the variation in orbit decay rate as a 
function of attitude control. Figure 1 7 shows the full 60,000-sec run duration, while figures 1 8 and 1 9 show 
shorter time segments taken near the end of the run. Figures 17-19 illustrate again the agreement of recur- 
sive and conventional solution methods for cases 1 and 2, while showing the method-induced deviation 
apparent in case 3, the no control (tumble) case. 

Appendix A.4 contains 42 plots that compare results from the two recursive approaches (RF-C and 
RF-D) to each other for all three cases. These plots are presented primarily for completeness since, for all 
practical purposes, the responses are identical for both methods over the full 60,000-sec duration of the 
simulation runs. 
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Figure 19. Altitude, case-method array, second time segment (km). 




8. PERFORMANCE 


The numerical results presented in figures 8-19 and in appendix A very clearly confirm that cou- 
pling of discipline-oriented subsimulations by recursive feedback can be used successfully as an architec- 
ture and solution approach to simulation of the example system. In cases I and 2, 60,000 sec of run time is 
achieved with no significant deviation of results. In case 3, deviation begins at approximately 45,000 sec. 
The earlier deviation of case 3 is attributed primarily to linear interpolation of aerodynamic forces and 
torques over the subsimulation integration step size in the recursive (RF— C and RF— D) implementations, 
while aerodynamic computations are made at points internal to the integration interval in the unified (RK- 
U) implementation. Improvement could likely be found through a more sophisticated interpolation method, 
a shorter convergence interval, or an increased number of points per convergence interval for signal defini- 
tion. It must be remembered, however, that eventual disagreement of dynamic simulations driven by differ- 
ent numerical processes is always expected. 

The results also confirm that the location and distribution of coupling functions, while partially 
defined by the fundamental nature of the subsystem interactions, can be part of the simulation design and 
can have a significant affect on performance. For the current example, the convergence interval length 
(0.4 sec) and the number of points per convergence interval (five) are the same for both recursive 
approaches, and the accompanying subsimulations as well as the conventional simulation all use the same 
integration step size (0.1 sec). Convergence error tolerances were in the range of 10“ 12 to 1() -14 for the 
recursive methods. Typical reported processor utilization times for the RF-C, and RF-D simulations are 
approximately 1,200 and 800 sec, respectively, while the RK-U time is approximately 200 sec. 

As mentioned earlier, the potential parallel computation advantage of the RF-C simulation is lost 
because of the particular single-computer FORTRAN implementation. Given that fact, the RF-D (fig. 6) 
simulation is understandably faster than the RF-C (fig. 4) because the amount of data being processed and 
checked to confirm or deny convergence at the recursive transfer point is significantly less. The speed of 
the RF-C simulation would be increased if the subsimulations were actually run in parallel. The possibility 
for running subsimulations in parallel does not exist for the RF— D simulation because the process design 
requires them to be run sequentially. In either recursive method, output from the last stage of the process is 
used to compare against the previous stage for confirmation of convergence. In the conventional method, 
no active control of error is currently present, and inclusion of any such mechanism would require addi- 
tional computation; conversely, predetermination of the number of recursion stages and elimination of the 
convergence check could significantly reduce the computation load for either recursive method. For the 
current example, the maximum number of recursion stages is typically three, one of which could arguably 
be said to serve no purpose other than to confirm satisfaction of error tolerances. As the number of recur- 
sion stages increases, the extra time required for the confirmation stage becomes less significant in a rela- 
tive sense, but elimination of the checking process entirely would save time at every stage. 
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None of the current simulations were constructed with emphasis on efficiency. After correctness, 
programmer convenience and listing readability were the next higher priority considerations. It is believed 
that many strengths of recursively coupled subsimulations architecture lie in areas other than computa- 
tional efficiency and those strengths may sometimes be well afforded at considerable expense. 
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9. IMMEDIATE FUTURE DEVELOPMENT 


More analytical investigation and development are appropriate. Analyses that investigate error 
sources, error propagation, computational efficiency, and stability are needed. Analyses that address appli- 
cability to stochastic systems and discrete data systems are also needed. Further reduction of conservatism 
ot the criterion for convergence derived for systems with algebraic loops (ref. 6) should be explored. The 
basic concept should be extended to include automatic control of convergence interval duration and other 
process parameters to improve performance. Implications of central versus distributed or mixed coupling 
architectures should be considered from viewpoints of algorithm performance and computational resource 
utilization. 

To facilitate a multicomputer implementation, a high-level standard format for digital definition of 
a vector analog signal segment (function of time) needs to be chosen or created so that signal segments can 
be transmitted and received throughout the system. In addition, system-level command and control func- 
tions must be defined so that standard subsimulation and coupler interfaces can be developed. A single- 
computer multiprocessor development that takes advantage of parallel computation could advance devel- 
opment of requirements for a multicomputer implementation as well as facilitate a better assessment of 
performance. Construction ol an example involving more than one computer will serve to identify a de- 
tailed set of platform- and network -specific requirements in data systems and software terms; it will also 
serve to explicitly demonstrate technical integrity for a multiplatform application. 
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10. CONCLUSIONS AND RECOMMENDATIONS 


An example nonlinear, continuous simulation of satellite motion has been successfully constructed 
using two variations of recursive feedback to couple three separate subsimulations of orbit dynamics, 
attitude dynamics and control, and aerodynamics. Results have been verified by direct comparison with a 
conventionally constructed simulation ot the same system. Each subsimulation deals with one engineei ing 
discipline, and appropriate interactions are implemented recursively at the system level. The subsimulations 
can, in principle, be run separately on remote and dissimilar platforms while coupling is implemented by 
way of network. Depending on the system-level design configuration, it is possible to run the subsimulations 
in parallel. The coupling signals represent physically identifiable variables associated with physical inter- 
faces of the subsystems. The example employs three subsimulations; the method framework accommo- 
dates n. Clearly, further investigation of recursively coupled subsimulations as a systems simulation 
architecture is warranted. It is recommended that the present effort be expanded to include all of the future 
development areas mentioned above, and that longer range plans should include a pilot multisubsystem, 
multidiscipline, multiplatform development of a complexity level that demonstrates effectiveness and 
efficiency through collaboration. 
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APPENDIX A— RESPONSE PLOTS 


A.l Response Plots for Case 1 (Full Three-Axis Inertial Attitude Control) 
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A.2 Response Plots for Case 2 (Two- Axis Attitude Control) 
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A .3 Response Plots for Case 3 (No Attitude Control) 
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A.4 Response Plots for Cases 1-3 With Central and Distributed Coupling 
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APPENDIX B— FORTRAN LISTINGS 


B.l FORTRAN Listings, Including Main Program and All Subroutines, 
Recursive Feedback With Centralized Coupling (RFC) 



This is FORTRAN main program "rfc.f". Development began 6/12/2000 by Doug Mixon, MSFC/TD55, 
in support of AEI project, and the listing is contained in "aei/rfc" 
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point for the vector output signal from the first (orbit dynamics) 
and is also the first (top) partition of vy(i,j). 
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print* 

print* , 'This is main program rfc.f. Solution 
print* , 1 method is recursive feedback with' 
print* centralized coupling.' 



print* , Input the number of data (time) poi 
print* , ' per convergence interval, including 
print*,' both end (boundary) points:’ 
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tol= (1 .ODO)* 1 .Oe-12 ! relative 

res= (1. 0 DO )* 1 . Oe-12 ’absolute 

itmax=100 Imaximum allowable number of iterations 

ptol=1.0e-3 ! tolerance for print time 

print* , 1 Run in progress.' 
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Decomposing the system-level feedback vector into individual subsystem input vectors 



£ 


a> 


P 


to 

i — t 

>1 

03 

60 

> 

X 

U 

3 

a 

to 

p 


c 

JH 

•rH 

a> 


a 

a 


o 

to 

c 

f — i 

a 

3 

co 

C 

M 

Co 

a 

•r-J 

> 

CO 

c 


o 

P 

u 

3 


a 

a 

c 

X 


p 

u 

p 

at 

a 

rH 

> 

to 

0 

u 


to 

CO 


c 

X 

o 

* H 

-rH 

co 

P 

— ■ 

05 


■ — t 

. — „ 

3 

O' 1 

E 

3 

■H 

3 

(0 

CM 

E 

3 

a 

3 

p 


CO 

rH 

!>i 

3 

CO 

3 

X 

«. 

3 

n 

to 

3 


V 

rH 

Qu 

03 

X 

3 

C 

X 

«, 

* rH 


> 

- — ' 

*r) 

a 

X 

E 

C 

o 

•H 

o 


at 

CO 

X 

c 


•rH 

< — 1 

c 

1 — 1 

c 

3 

3 

U 

Qd 


* 

o o o u 


* 

¥ 

¥ 

* 




03 

* 






> 

¥ 






U 

¥ 






a 

¥ 






p 

¥ 






c 

¥ 





t — 1 

-H 

¥ 





03 


¥ 





> 

a 

¥ 




X 

M 

0 

¥ 




03 

a 

c 

¥ 




> 

p 

a 

¥ 




M 

c 

Co 

¥ 




a 

•H 

u 

¥ 


Co 


p 


a 



G 


c 

a 

> 

¥ 


■rH 


■ w 

u 

c 

¥ 


1 — 1 



a 

0 

¥ 


a 


a 

a 

0 

¥ 


3 


o 

co 


¥ 


0 


c 

u 

u 

¥ 


0 


a 

a 

a 

¥ 




Do 

> 

> 

¥ 


E 

CO 

H 

c 

0 

¥ 


a 

a 

a 

o 


¥ 


p 

3 

> 

u 

c 

¥ 


CO 

rH 

c 


0 

¥ 


>1 

03 

o 

Jh 

•H 

¥ 


to 

> 

o 

a 

p 

¥ 


X 



> 

03 

¥ 


3 

p 

u 

o 

rH 

* 


(0 

C 

a 


3 

¥ 



a 

> 

c 

E 

¥ 



u 

o 

o 

-H 



1 — 1 

Jh 


■H 

CO 

P 


03 

3 

c 

p 


3 


G 

U 

o 

03 

. — . 

a 


Co 


•H 

i — 1 

a 

p 


• rH 

O 

p 

3 

•r-f 

3 


co 

P 

03 

E 

03 

0 




<X 

■rH 

Jh 



X 

03 

3 

CO 

X 

rH 


O 

P 

E 


a 

a 


03 

03 

•r-f 

CO 

co 

> 


X 

X 

CO 

u 

rH 

a 


X 



•H 

03 

rH 


a 

a 

CO 

E 

■ — ' 

1 


a 

Co 

u 

05 


E 


p 

03 

■H 

C 

CO 

a 



p 

E 

£0 

u 

p 


rH 

CO 

03 

X 

•rH 

CO 


a 


C 


E 

>1 


> 

c 

>i 

a 

03 

CO 


a 

0 

T 3 

X 

C 



p — 1 

■rH 


3 

>1 

a 


1 

CO 

P 

p 

X 

1 — 1 


E 

Jh 

-«H 

•rH 

0 

co 


a 

3 

X 

P 

u 

c 


p 

O 

JH 

P 

a 

■H 

>, 

CO 

a 

0 

03 

03 

CO 

>1 

>t 

jh 






CO 


C 

C 

c 

O 

no 


1 — I 

3 

3 

3 

P 

>1 

CO 

3 

M 

M 

M 

c 


a 

■H 




*H 

*. 

p 

P 


- 

•w 


CM 

3 

■H 

on 

00 

X 

CO 

>1 

a 

C 




P 

>1 

E 

■rH 

a 

a 

a 

3 


0 


X 

X 

X 

a 

iH 

0 

CO 

o 

0 

0 

p 

rO 


p 

E 

E 

6 

3 

>1 


a 




0 

h. 


CO 

— • 

— • 

— • 


a 

. — 

— . 




E 

X 

3 





a 

c 

3 



- — - 


p 



_ _ 


no 


CO 

<D 

>1 

* 1 — \ 


' — ' 


>1 

— ' 

>1 



a 

— 

CO 

a 

- — - 

•H 

00 

X 

c 

X 

H 

u 

a — 

— ' 

3 

>1 

3 

X 

a 

00 X >1 

P 

p 

X 

CO 

E 

rH 

H G >1 

*H 

*H 

0 


a 

a 

^ 11 

X 

P 

u 

a 

CO 

3 

t 1 ' 1 —V d) 

sn 

P 

a 

* — I 

CO 

0 

II II ■( t 3 

o 

03 

03 

X 

03 

0 

•H -n C 




E 



■H -H 

l-H 

< — 1 

i— t 

a 

1 — 1 

> — 1 

2 

2 

( 

t 

1 — t 

• — I 

rH 

CO 

• — 1 

1 — t 

a c 

03 

03 

03 

CO 

03 

03 

0 0 >1 0 

O 

O 

O 

< 

U 

O 

X X >1 O 




¥ 



CM 




-K 

u U U CJ u u 


+ * * 

* * 4c 

* * ¥ 

* * * 

¥ ¥ ¥ 

+ + * 

¥ -K * 


¥ 

¥ ¥ 






1 — 1 

¥ 

¥ ¥ 






3 

¥ 

¥ ¥ 






j> 

¥ 

¥ ¥ 






Jh 

¥ 

¥ ¥ 






a 

¥ 

¥ ¥ 






p 

¥ 

¥ ¥ 






c 

¥ 

¥ ¥ 





X 

■H 

¥ 

¥ ¥ 





3 


¥ 

¥ ¥ 

CO 




> 

a 

¥ 

¥ ¥ 

u 



rH 

Jh 

0 

¥ 

¥ * 

0 




a 

c 

¥ 

+ ¥ 

p 



> 

p 

a 

¥ 

¥ ¥ 

0 



Sh 

c 

co 

¥ 

¥ ¥ 

a 



a 


Jh 

¥ 

¥ ¥ 

> 



p 


a 

¥ 

¥ ¥ 




G 

a 

> 

¥ 

¥ ¥ 

p 



•H 

0 

G 

¥ 

¥ ¥ 

3 




c 

O 

¥ 

¥ ¥ 

a 



a 

a 

O 

¥ 

¥ ¥ 

G 

E 


V 

Co 


¥ 

¥ ¥ 


a 


c 

Jh 

Jh 

¥ 

¥ ¥ 


p 


a 

a 

a 

¥ 

¥ ¥ 

£ 

CO 

1 — 1 

Co 

> 

> 

¥ 

¥ ¥ 

a 

D >1 


u 

c 

0 

¥ 

¥ ¥ 

p 

CO 

> 

a 

0 


¥ 

¥ ¥ 

CO 

X 

M 

> 

u 

G 

¥ 

¥ ¥ 


3 

a 

c 


0 

¥ 

¥ ¥ 

CO 

CO 

p 

0 

M 

'H 

¥ 

¥ ¥ 

X 


G 

0 

a 

P 

¥ 

¥ ¥ 

3 


■rH 


> 

3 

¥ 

¥ ¥ 

(0 

a 


in 

0 

rH 

¥ 

¥ ¥ 


a, 

a 

a 


3 

¥ 

¥ * 

1 — 1 


0 

> 

c 

E 

¥ 

¥ ¥ 

3 

CO 

c 

0 

0 

•H 

¥ 

¥ ¥ 

3 

1 — 1 

a 


*rH 

CO 

¥ 

EH ¥ 

X 

3 

co 

G 

P 


¥ 

X ¥ 

♦H 

G 

Jh 

O 

3 


¥ 

X * 

> 

O' 

a 

•rH 

X 

0 

¥ 

O * 

•H 

•n 

> 

P 

3 


¥ 

CU + 

X 

CO 

c 

3 

E 

3 

¥ 

¥ 

G 


0 

r-H 

-H 

Jh 

¥ 

^ * 

•r-i 

P 

0 

3 

CO 

X 

¥ 

PC ¥ 


3 


E 


a 

¥ 

Eh * 

O 

a 

a 

•rH 

CO 

CO 

¥ 

^ + 

P 

G 

X 

CO 

O 

X 

¥ 

W + 

c 

•rH 

p 


•H 

3 

¥ 

¥ 

•H 



CO 

E 


¥ 

P-i ¥ 


5h 

VH 

u 

3 


¥ 

O + 

M 

a 

a 

•H 

C 

CO 

¥ 

O ¥ 

O 

rH 

> 

E 

>1 

O 

¥ 

X He 

P 

3 

0 

3 

X 

■rH 

¥ 

¥ 

U 

O 


G 


£ 

¥ 

* 

a 

CO 

CO 

>1 

a 

3 

¥ 

O * 

> 


c 

X 

X 

G 

¥ 

H ^ 


K 

0 


3 

>1 

¥ 

CO ¥ 

M 

•H 

■H 

P 

p 

X 

¥ 

PC ¥ 

O 

(0 

P 

•H 

•H 

O 

¥ 

O ¥ 

3 

— ■ 

3 

X 

p 

Jh 

¥ 

O ¥ 

X 


rH 

Jh 

p 

a 

¥ 

W hc 

X 


3 

O 

3 

3 

¥ 

PC ¥ 

a 

on 

E 




¥ 

¥ ¥ 

a 

3 

-H 

C 

G 

G 

¥ 

¥ ¥ 

p 

3 

CO 

3 

3 

3 

¥ 

¥ ¥ 




Jh 

Jh 

Jh 

¥ 

¥ ¥ 

H 

CM 

E 




¥ 

¥ ¥ 

a 

3 

a 




¥ 

¥ ¥ 

> 

3 

p 

on 

on 

on 

¥ 

¥ ¥ 

a 


CO 




¥ 

¥ * 

1 — t 

rH 

»>i 

a 

a 

a 

¥ 

* ¥ 

1 

3 

CO 

X 

X 

X 

¥ 

¥ ¥ 

E 

3 

X 

0 

0 

0 

¥ 

* ¥ 

a 


3 

E 

E 

E 

¥ 

¥ ¥ 

p 

3 

co 




¥ 

¥ ¥ 

CO 

3 



_ . 


¥ 

¥ ¥ 

>1 

N. 

1 — 1 




¥ 

¥ ¥ 

CO 

a 

3 




¥ 

¥ ¥ 


X 

3 




¥ 

¥ ¥ 

a 

c 

X 


on 


¥ 

¥ ¥ 

X 

•s 

-H 


- — - 

on 

¥ 

¥ ¥ 

p 


> 


a 

— 

¥ 

¥ ¥ 



•rH 

on 

X 

G 

¥ 

¥ ¥ 

at 

a 

X 

- — 

3 

>1 

¥ 

— • He 

c 

£ 

c 

P 

p 

X 

¥ 

¥ 

•rH 

0 

•H 

■rH 

*H 

O 

¥ 

a He 

CO 

u 


X 

p 

U 

¥ 

3 He 

0 

a 

Co 

Jh 

p 

a 

¥ 

C ¥ 

a 

X 

c 

O 

3 

3 

¥ 

•H He 

E 


-H 




¥ 

P He 

0 

< — f 

c 

1 — 1 

X 

• — 1 

¥ 

C He 

0 

f — 1 

c 

1 — 1 

1 — 1 

X 

¥ 

O He 

a 

3 

3 

3 

3 

3 

¥ 

O * 

Q 

0 

PC 

U 

O 

O 

¥ 

He 







¥ 

O He 

¥ 


¥ 




¥ 

O He 

¥ 


¥ 




¥ 

CM He 

¥ 


¥ 




O 

c 

c 

O O 

u 

C 

C 

O 

O 



96 


Assemble subsystem outputs into single system-level output 
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dtmaj ! simulated problem time , for reference only 
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if (nmt . eq . nmtr ) print*, ’ Normsl 0xit, run-time = 1 , irun, * sec 
print* 
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dimension cpf ront ( 3) , cpright ( 3) , cpbottom (3) 
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dimension a front (3) , aright (3) , abottom (3 ) 



dimension aproj (6) , drag (6) , face for (6,3), facet or (6, 3) 
dimension t force ( 3 ) , ttorq ( 3 ) 
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facefor(i,j) drag ( i ) * vnorm^b ( j ) ! this is jth component of force on the ith 



22 continue ! facet, expressed in body coordinates 
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do 36 i=l, 3 ! adding rear facet contribution 

t force (i)=facefor(4, i) 


call vcross (cprear , tf orce, ttorq) 

do 37 i=l , 3 

facetor ( 4 , i) =ttorq(i) 
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subroutine aerodyn (mode ) ! mode 1= set physical constants, system parameters, etc. 

! mode 2= blank, no initial conditions for algebraic system 
! mode 3= run simulation over the convergence interval 
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100 continue ! (mode-1, define constants, system parameters) 
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af ront ( 1 ) =ax ! setting non-O.ODO components of facet area vectors 
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End of mode 2 , beginning of 



300 continue ! (models, run simulation over the convergence interval) 
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uu2 (6, 100) , yy2 (6, 100) ! I /O for subsystem 2 (attitude dynamics) 

yy2i (6) , yy2f (6) !; initial and final output conditions, subsystem 
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Gravity constant (gme) is defined by orbit. f, and is obtained through bloc 
common /constants/; pi (3.14xx...) and cvdr (deg. to rad. mult, factor) an 
are defined and obtained similarly. 
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print* , ' (integer) : 

read* , idiv 
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phi (i) =yy2i (i) 
omega (i) =yy2i (i+3) 

yy2 (i, 1) =phi ( i) (entering initial output vector (time point number 1) into common/subsys2 / 

y y 2 ( i + 3 , 1 ) =omega ( i ) 
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f yy / yyp / i conf 1 ag ) ! system level convergence check 
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from orbit dynamics to attitude dynamics: 



rad ( 1 ) =yy ( 1 , j ) [Extracting orbit radius from the system output signal 

rad (2 ) =yy (2 , j ) [variable rad(i) is used for listing readability only 

rad (3) =yy (3, j ) 
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vel ( 1 ) =yy ( 4 , j ) [variable vel(i) is used for listing readability only 


call matvec (amat, vel, vel_b) [transforms orbit velocity vector to body coordinates 
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phi (3) =datan2 (amat (1, 2) , amat (1, 1) ) ! (+/- p i 
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continuing search for the interval 



performing the interpolation 


Q) 

1 * o 

c 

O -H 

4-> 4-J 

C 

o o 

tn o 

CO 


o 





X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 



>— * 

X 

X 

X 

X 

X 

X 



•H 

X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 



o 

X 

X 

X 

X 

X 

X 



c 

X 

X 

X 

X 

X 

X 



D 

X 

X 

X 

X 

X 

X 



4-1 

X 

X 

X 

X 

X 

X 



f 

X 

X 

X 

X 

X 

X 



. — . 

X 

X 

X 

X 

X 

X 



i — I 

X 

X 

X 

X 

X 

X 



+ 

X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 



' — ' 

X 

X 

X 

X 

X 

X 



u 

X 

X 

X 

X 

X 

X 



c 

X 

X 

X 

X 

X 

X 



D 

X 

X 

X 

X 

X 

X 


-l-> 

4-1 

X 

X 

X 

X 

X 

X 


i — 1 

• — ' 

X 

X 

X 

X 

X 

X 


Q) 

4c 

X 

X 

X 

X 

X 

X 


V 

o 

X 

X 

X 

X 

X 

X 


\ 

ro 

X 

X 

X 

X 

X 

X 



4-1 

X 

X 

X 

X 

X 

X 


i — i 

+ 

X 

X 

X 

X 

X 

X 


4-> 


X 

X 

X 

X 

X 

X 

» — 1 

1 

■H 

X 

X 

X 

X 

X 

X 


0) 

' — ' 

X 

X 

X 

X 

X 

X 

> 

£ 

O 

X 

X 

X 

X 

X 

X 

U 


c 

X 

X 

X 

X 

X 

X 

Q) 

-M 

p 

C X 

X 

X 

X 

X 

X 

4-> 

■ — ■ 

4-j 

M X 

X 

X 

X 

X 

X 

C 

II 

tl 

P X 

X 

X 

X 

X 

X 

■H 

U 

i — 1 

■P TJ X 

X 

X 

X 

X 

X 

1! 

m 

rO 

(DC X 

X 

X 

X 

X 

X 

-H 

mh 

> 

MO) X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 




X 

X 

X 

X 

X 

X 



O 

X 

X 

X 

X 

X 

X 



X 















M 















pc 















H 















< 















£ 















ro 



s 












X 



ro 












ro 



hh 



« — .. 

. — . 

.—l 

.. — i 


, — i 


, — . 





CM 

— V 


ro 

ro 

CM 

CO 

CO 

04 

CO 

CO 


< 




ro 

CO 


hi 

hi 

1 








- — - 

< 

hi 


CM 

CO 

CM 

CO 

rH 

CO 

! 1 

04 


Pu 

CO 

ro 

4C 

T 1 

1 — 1 

— " 

- — ' 

- — 

■ — 

' — 

' — ^ 





o 

1 

■h 


•— ' 


< 

< 

< 

< 

< 

< 

< 

< 



o 

ro 

CM 

< 

< 

4C 

4c 

4c 

4c 

4c 

4c 

4c 

4c 


w 

•h 

' — ' 


4C 

4< 


. — l 



,_i 

^ i 




CO 

X! 

Pm 

ro 



CM 

1 — 1 

1 — 1 

CM 

rH 

rH 

04 

i — 1 


PC 

1 

£ 


CM 

CM 

h. 

h. 

h. 

ii 

h. 

,i 


H 


u 

< 

pa 

< 

h. 

hi 

ro 

CM 

CO 

rH 

CO 

rH 

04 

T 1 


> 


Eh 

t 

ro 

C4 

— 

— 

- — ' 

— ' 

— - 

> — ' 

•• — - 



z 


V 




< 

< 

< 

< 

C 

< 

< 

< 


M 

o 
1 — 1 

ro 

ro 

<c 

< 

1 

1 

1 

1 

i 

1 

1 

1 


tx] 

CO 


ro 


, 

ro 

ro 

CM 

00 

CO 

04 

CO 

ro 

—v 

DC 

M 

ro 

" — " 

ro 

ro 


h. 


hi 

w 




M 

H 

u 

- — 

< 

hi 

hi 

ro 

CM 

ro 

I — 1 

CO 

rH 

04 

x — 1 

<c 


pa 

M 

4C 

ro 

CM 

' — ' 

— ' 

- — ' 


- — - 


— - 


s. 

CO 

PC 

< 




< 

< 

< 

< 

< 

< 

< 

< 

< 

w 

CM 

•h. 

CM 

<c 

< 

4c 

4c 

4c 

4c 

4c 

4c 

4e 

4c 

' — 

X 



h. 

4c 

4c 


, — l 

i— i 






CO 

< 

pa 

ro 

CM 



04 

1 — 1 

rH 

CM 

rH 

rH 

04 

rH 

> 

Eh 

bC 

V 

— 

CM 

04 


h. 

h. 

h. 

h. 




£ 


PQ 

ro 

< 

V 


CM 

CO 

CM 

CO 

rH 

CO 

i — 1 

04 

M 

w 

ID 

' — ' 

■ — 

1 — 1 

rH 

' — ' 

' — 

■ — 

' — 



— - 

— 

— - 


£ 

Q 

< 

* 



< 

< 

< 

< 

< 

< 

< 

< 

w 

1 — 1 

Q 


. — l 

< 

< 

II 

II 

II 

Ii 

II 

II 

II 

II 


Eh 



rH 

— 

' — ' 



, — i 

, — l 



, i 

„ i 

w 

ID 

Eh 

o 

*h 

4C 

4C 

rH 

CM 

CO 

< — 1 

CM 

CO 

rH 

04 

H 

9 

1 — 1 

M 

rH 



hi 

>h 

V 

V 



v 


D 

PC 

u 

CO 

— ' 

< — 1 

rH 

1 — ( 

rH 

rH 

04 

04 

04 

CO 

00 

o 


M 

2 

< 

n 



' — - 

— " 

' — ■ 



' — 

' — 

pc 

CO 


pa 

ii 

04 

ro 

cu 

CM 

a. 

CM 

Cm 

04 

Pu 

PM 

PQ 

H 

Pm 

2 

EH 



£ 

2 

£ 

S 

£ 

2 

£ 

£ 

D 

a: 

£ 

M 

pa 

< 

< 

ca 

pa 

pa 

pa 

pa 

pa 

pa 

pa 

CO 

Fh 

M 

a 

Q 

1 

+ 

H 

Eh 


Eh 

Eh 

H 

Eh 

Eh 


o o 


4c 

-K 

4c 

u u u u o 


122 


X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X x 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X x 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 
X X X X X X 







X 

X 

X 

X 

X 

X 


X 














X 

X 

X 

X 

X 

X 


M 














X 

X 

X 

X 

X 

X 


PC 














X 

X 

X 

X 

X 

X 


H 














X 

X 

X 

X 

X 

X 


< 














X 

X 

X 

X 

X 

X 


2 














X 

X 

X 

X 

X 

X 











. 





X 

X 

X 

X 

X 

X 


< 









CN 





X 

X 

X 

X 

X 

X 
















X 

X 

X 

X 

X 

X 


X 









rH 





X 

X 

X 

X 

X 

X 


PC 

-—V 








— - 





X 

X 

X 

X 

X 

X 


PC 

CO 








< 





X 

X 

X 

X 

X 

X 


1 













X 

X 

X 

X 

X 

X 


o 

o 













X 

X 

X 

X 

X 

X 


H 









H 





X 

X 

X 

X 

X 

X 


CJ 

PC 













X 

X 

X 

X 

X 

X 


to 

1 








CNJ 





X 

X 

X 

X 

X 

X 


> 

< 













X 

X 

X 

X 

X 

X 



" — ■ 

00 







< 





X 

X 

X 

X 

X 

X 


< 

2 








1 





X 

X 

X 

X 

X 

X 

> 


o 

> 












X 

X 

X 

X 

X 

X 

< 

CO 

M 

< 







CsJ 





X 

X 

X 

X 

X 

X 

< 

to 

CO 













X 

X 

X 

X 

X 

X 

> 

M 

M 

- — - 







CN 


CO 



X 

X 

X 

X 

X 

X 



u 

00 




> 





Q 



X 

X 

X 

X 

X 

X 

< 

CM 

to 

' — " 












X 

X 

X 

X 

X 

X 

' — ' 

\ — \ 

PC 

> 




' — ' 



* 






X 

X 

X 

X 

X 

X 

(J 

H — 

On 





•d 








X 

X 

X 

X 

X 

X 

to 

►0 T— 1 






*■“ 



< — 1 





X 

X 

X 

X 

X 

X 

> 

D X 

to 

00 




1— 1 








X 

X 

X 

X 

X 

X 


2 00 

►0 








T \ 





X 

X 

X 

X 

X 

X 

< 

' — " 

CQ 

00 




< 





Cm 



X 

X 

X 

X 

X 

X 

s 

tO II 

D 

' — ' 




+ 



< 

OO ro 

2 



X 

X 

X 

X 

X 

X 


2 ~ 

Q 

< 

00 


oo 




II 


to 



X 

X 

X 

X 

X 

X 

to 

1 — 1 »— 1 

Q 





1— 1 




r-H i — 1 

E-h 



X 

X 

X 

X 

X 

X 

2 

H X 


2 

r- H 

o 

1 — 1 

' 



ro 

II l! 

II 



X 

X 

X 

X 

X 

X 

1 — 1 

D oo 

E-« 

o 

II 

• 

'1 

> 




M 




X 

X 

X 

X 

X 

X 

H 

o ^ 

H 

1 — t 

M 

o 


< 



oo 


»“> 

2 


X 

X 

X 

X 

X 

X 

D 

PC + 

o 

CO 


II 


II 




O O 


PC 


X 

X 

X 

X 

X 

X 

o 

— 

1 — i 

2 

o 

- — - 

o 

' — - 

PC 


Cm 

< — 1 i — 1 

M 

D 


X 

X 

X 

X 

X 

X 

PC 

CO 00 

kO 

to 

1 — 1 

M 

rH 

M 

D 


2 


- — ' 


Q 

X 

X 

X 

X 

X 

X 

CQ 

M X 

On 

2 


' ' 



E-I 

Q 

to 

H 

88 

M 

c 

to 

cc 

to 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

D 

CO 

X 00 
H ~ 

2 

M 

M 

Q 

DO 

> 

< 

8 

> 

< 

to 

PC 

2 

to 






X 

X 

X 

X 

X 

X 













o 



X 

X 

X 

X 

X 

X 


* 






o 





\ — i 



X 

X 

X 

X 

X 

X 


* 






T 1 








X 

X 

X 

X 

X 

X 


* 














X 

X 

X 

X 

X 

X 





V 











X 

X 

X 

X 

X 

X 


c 

c 

c 

u 

u 







X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


123 



subroutine orbit (mode) ! mode 1= set physical constants, system parameters , etc 

! mode 2= set initial conditions 

! mode 3= run simulation over the convergence interval 
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Numerically Defining Constants (mod 



Constant of Gravitation, 
m* * 3 /kg- sec** 2, CRC, page 
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200 if (mode. eq. 3) go to 300 

continue ! (mode = 2, set initial conditions) 



Setting Initial Conditions (mode 2 function) : 
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300 continue ! (mode-3, run simulation over the convergence interval) 


! r i 0 v ing initial conditions from block common /endcondl/ 
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end integration loop 
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!post processing orbit data for altitude 



call vecnorm (rad, rad_n, rmag) 

alt=rmag/1000 . ODO-rearth laltitude, kilometers 
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IMPLICIT DOUBLE PRECISION (A- 



THIS SUBROUTINE PERFORMS ONE STEP IN THE RUNGE-KUTTA 
INTEGRATION OF A SYSTEM OF FIRST ORDER DIFFERENTIAL 
EQUATIONS OF THE FORM Y ' (X) =FD (X, Y) , WHERE Y AND FD 
ARE VECTOR FUNCTIONS, AND X IS THE INDEPENDENT 
VARIABLE. THE ALGORITHM IS DESCRIBED BY THE FOLLOWING 
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CALL SUB1 (FDERIV, Y, X) ! "1" added to identify with subsystem #1 

DO 500 1=1, N 



YTEMP ( I ) =YTEMP ( I ) +FDERIV ( I ) * H3 
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INTEGRATION OF A SYSTEM OF FIRST ORDER DIFFERENTIAL 
EQUATIONS OF THE FORM Y ' (X) =FD (X, Y) , WHERE Y AND FD 
ARE VECTOR FUNCTIONS, AND X IS THE INDEPENDENT 
VARIABLE. THE ALGORITHM IS DESCRIBED BY THE FOLLOWING 
VECTOR RELATIONS: 



F1 = FD (X, Y) , YT=Y + F1* H/6 

F2="FD(X+H/2,Y + Fl*H/2) , YT=YT+F2*H/3 
F3=FD(X+H/2,Y+F2*H/2) , YT=YT+F3*H/3 
F4=FD (X+H, Y+F3* H) , ' Y=YT+F4*H/6 
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X=X0+H 

CALL SUB2 (FDERIV, Y, X) ! "2" added to identify with subsystem 

DO 600 1=1, N 
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xdot (2 ) -rdot (2 ) 
xdot (3) =rdot (3) 
xdot (4) =vdot (1) 
xdot (5 ) =vdot (2) 
xdot ( 6 ) =vdot ( 3 ) 
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total ( 3 ) ! total external torque minus gyroscopic torque 
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x(l)=phi(l): Euler angle rotation from inertial to body system, about x-axis, 
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call vecnorm ( rad, rad_norm, rmag) ! normalized radius vector in inertial coordinates 
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control (2) =- (bO* phi (2) +bl* omega (2 ) ) 
control (3) =- (cO* phi (3) +cl* omega ( 3 ) ) 
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Changing back to mathematical notation from physical notation. 



xdot (2 ) ^phidot ( 2 ) 
xdot (3 ) =phidot ( 3 ) 
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Third subsystem output signal. Represents the ith component at the jth time 
point for the vector output signal to the third (aerodynamics) subsystem, and is 
also the third (bottom) partition of yy ( i , j ) . 



yyli(i) Value of the output signal at the beginning of the 
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of the distributed coupling configuration, information flows from one sub-simulation to the 
a sequential manner. The order of execution is important, and parallel execution is not 
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. 0/ float (ndp-1 ) ! ratio of integration time step to convergence interval duration 



print* , ' The above data limits all integration' 
print*,' stepsizes to ',delt,' sec.' 
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print*,' Run in progress.’ 
print* 

call clock (nl ) Imachine start time 
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Running individual subsystem simulations over the convergence interval 

call orbit (3) ! mode 3, run orbit dynamics simulation over convergence interval 


call attitude (3) ! mode 3, run attitude dynamics simulation over convergence interval 





rH 


(J 

ID 

CO 

iD 

> 

o 

XI 

P 

•«H 

73 

a; 

e 

0 ) 

p 

fD 

at 

c 

c 

p 

■H 

>1 



X 3 


a> 


CO 

u 

a> 

o 

c 

T 3 

•H 

0 ) 

3 

£ 

CP 

P 

<D 

p 

■H 

G 

at 

P 

>i 

> 

p 

T 3 

c 

iD 

O 

o 


P 

o 

73 

at 


C 

ID 


P 

0 ) 

> 

O 

c 

o 

-H 

P 

ID 

<-1 

3 

6 


r0 

P 

XI 

Qt 

CP 


XI 

p 

o 

o 

P 

CO 

u 


ID 

c 

T3 

O 

P 

Q) 

rO 


CO 

0 t 

3 

rH 

1X3 

> 

P 

c 

at 

x 

P 

3 

o 

o 

P 

1X3 

p 

ID 

"O 


(0 

a> 

P 

1 X 3 CO 
C (U 

-H p 

73 ID 
5-1 C 
O *H 
O X3 
a p 
o 

> — ( o 

iD U 


ID 

O 

a 

p 



P 

CP 

p 

X 3 


P 

ID 

at 

o 

CO 


P 

c 

X 3 

a 

rH 

CO 

-rH 


■H 

ID 




£ 

C 

c 


a> 

fD 

CP 

o 

at 

3 

C 

“ rH 


o 

CD 


CO 

(0 

p 

P 

*o 


p 

o 

o 

o 


3 

p 

p 

p 

u 

(J 



a) 

tD 

at 

u 

o 

rD 

XI 

p 

•H 

-H 


T 3 


£ 

£ 

G 

at 

• — i 

<D 

<D 

3 

at 

ID 

G 

C 

P 

p 

-H 

>i 

>i 



P 

T 3 

V 


CO 

•rH 

O 

o 

00 

at 

c 

P 

p 


p 


at 

0 ) 

at 

3 


ID 

iD 

T 3 

a 

CO 



O 

£ 

p 



£ 

O 

at 




u 

CO 


■n 


CM 

oo 

00 


m H 

P 

C 

p 


CP ' — i 

at 

>i 

a> 


T 3 3 

i— H 

T 3 

rH 

00 

C 3 

a 

O 

a 


- II 

3 

P 

3 

1 — 1 

iH * — .. 

O 

at 

O 

II 

II -m 

U 

iD 

O 

*H 

■n 

rH 

i — i 

rH 

CM 

•H 

CM 

' — 1 

rH 

*H 


a 

ID 

(D 

ID 

0 

O XI 

U 

o 

U 

X 3 

73 p 


OO 

+ 


CM 

3 

3 

II 


X 

x 

* 

* 

X 

X 

X 

X 

X 

■X 

* 

* 

X 

-X 

4 c 

•X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


CO 

0) 

CP 

ID 

p 

CO 

c 

o 

■H 

CO 

p 

3 

o 

0) 

p 

4-1 


ID 

> 

p 

CD 

p 

c 


ID 

> 

p 

at 

p) 

g 

-H 

at 

u 

c 

at 

CP 

P 

(D 

> 

c 

o 

o 


CO 

od 

D 

O 

W 

a: 


a) 

p 


0) 

x: 


p 

(D 


3 

E 


p 

-H 

JQ 

P 

O 


iD 

> 

P 

0) 

P 

C 

•H 

a) 

u 

c 

a> 

CP 

p 

a) 

> 

c 

o 

u 

p 

CD 

> 

o 


3 

P 

•H 

p 


ID 

> 

P 

CD 

P 

C 


CD 

o 

C 

Gt 

CP 

P 

at 

> 

G 

O 

o 

p 

a; 

> 

o 

c 

o 


X 

X 

X 

o 

0) 

at 


p 

X 

X 

X 


u 

> 

C 

E 

X 

X 

X 

p 

c 

o 

o 

■H 

X 

X 

X 

at 

at 


•H 

CO 

X 

E-t 

X 

£1 

CP 

c 

p 


X 

2: 

X 

£ 

p 

o 

iD 


X 

1 — 1 

X 

3 

at 

-H 

t — i 

o 

X 

o 

X 

C 

> 

P 

3 

* H 

X 

Q-t 

X 


c 

3 

E 

(D 

X 


X 


o 

rH 

■H 

u 

X 


X 

P 

o 

3 

(0 

£l 

X 

cc 

X 

G 


£ 


at 

X 

H 

X 

3 

at 

■H 

CO 

Co 

X 

2 

X 

O 

x 

CO 

u 

, — i 

X 

IX 

X 

u 

p 


■H 

ID 

X 


X 



CO 

£ 


X 

Cu 

X 

CO 

p 

o 

fD 


X 

O 

X 

c 

at 

•r"* 

C 

CO 

X 

o 

X 

O 

> 

£ 

>i 

(J 

X 

(J 

X 

■H 

o 

<D 

T 3 

*H 

X 


X 

P 


C 


£ 

X 

s 

X 

ID 

CO 

>i 

at 

fD 

X 

o 

X 

P 

C 

73 

T) 

C 


>1 

73 

O 

P 

a) 

fD 


X 

X 

X 


-H 

C 


C 


c 


X 

X 

X 

CP 

CO 

3 


3 


p 


X 

X 

X 

c 


P 


P 


p 


X 

X 

X 

•H 

£ 







X 

X 

X 

p 

a; 

< 






X 

X 

X 

id 

p 

oo 


oo 


oo 


X 

X 

X 

73 

CO 







X 

X 

X 

a 

> 

at 


at 


at 


X 

X 

X 

3 

CO 



73 


73 


X 

X 

X 

— ■ 

XI 

o 


O 


o 


X 

X 

X 


3 

£ 


£ 


£ 


X 

X 

X 


co 







X 

X 

X 



— * 




__ . 


X 

X 

X 


rH 







X 

X 

X 


CD 







X 

X 

X 


3 



, ^ 




X 

X 

X 


T 3 



00 


___ 


X 

X 

X 


•p 



' 


00 


X 

X 

X 


> 


t — 1 

at 

CM 


00 

X 

X 

X 


“H 

oo 

P 

73 

P 

c 

p 

X 

X 

X 


TS 


at 

3 

at 

>1 

at 

X 

— ■ 

X 


c 

p 

i — i 

P 

• — i 

73 

r— | 

X 


X 

< — 1 

•H 

■H 

a 

*H 

a 

0 

Cu 

X 

at 

X 

+ 


XI 

3 

P 

3 

p 

p 

X 

3 

X 

p 

CP 

p 

o 

P 

o 

at 

0 

X 

c 

X 

*H 

c 

0 

U 

fD 

U 

fD 

u 

X 

■H 

X 

c 








X 

p 

X 

II 

c 

> — 1 

i — 1 

i — ( 

1 — 1 

r — | 


X 

c 

X 

p 

C 

rH 

• — f 

rH 

1 — 1 

fH 

r"H 

X 

o 

X 

•H 

3 

ID 

iD 

(D 

(D 

^3 


X 

u 

X 

G 

QC 

O 

U 

U 

U 

U 

u 

X 


X 









X 

o 

X 


x 







X 

o 

X 


X 







X 

CM 

X 


X 








Q) 

O 

G 

a> 

CP 

p 

0 ) 

> 

c 

o 

u 

p 

o 

p 

o 

a> 

x 

o 


o 

TJ 


a u 


u u o 


u u u 


154 



force, inertial coordinates 
torque, body coordinates 
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if (nmt.eq.nmtr) go to 300 lend of run, time is greater than or equal to the stop time (t_stop) 
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END OF MAIN PROGRAM 
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uu2 (3, j ) =yyl (3, j ) 
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Coupling from orbit and attitude dynamics to 


rad ( 1 ) =yyl ( 1 , j ) lextracting orbit radius from subsystem 1 output 

rad (2) =yyl (2, j ) 
rad (3) =yyl (3, j ) 
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implicit double precision (a- 



dimension uul ( 6, 100) , yyl ( 6, 100) 
dimension uu2 (6, 100) , yy2 (6, 100) 
dimension uu3 (6, 100) , yy3 (6, 100) 
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zz(i+12,j)=yy3(i,j) 



phi (i) =zz (i + 6, j ) ! unlimited angle (10/2/00) 

call dcm321 (phi, amat) 
call dcm2ang321 (amat, phi) 

zz (i + 6, j)=phi (i) /cvdr ! limited angle (10/2/00) 

zz ( i+9 , j ) =yy2 (i+3, j) /cvdr 
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100 format (20E20 . 10, i5) 
time=t ime+delt. 
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B.3 FORTRAN Listings, Including Main Program and All Subroutines Not Previously 
Listed in A. I, Unified Approach, Runge-Kutta Method (RK-U) 
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print *, 1 This main program rku.f. It uses a’ 

print*,' conventional centralized integration' 

print* , ' method. (Runge-Kutta) . ' 

print* 

print* , ' Input the integration stepsize , 1 
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xmoi ( 1 , 2 ) =0 . 0D0 ! kilogram-meters* * 2 

xmoi ( 1 , 3) =0 . 0D0 
xmoi (2 , 1 ) =0 . 0D0 

xmoi (2 , 2 ) =1 93125 . 0D0 ! 3m diameter x 15m length uniform- 



xmoi (2, 3) -0 . 0D0 ! density right-circular cylinder 

xmoi ( 3 , 1 ) =0 . 0D0 
xmoi ( 3, 2 ) =0 . 0D0 

xmoi (3, 3) =193125. 0D0* .9 ! IZZ=. 9* Iyy 
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cpf ront ( i ) — 0 . 0D0 ! centBrs of pressure location vectors ; from geometric 

cprear (i) -0 . 0D0 ! center of box to center of each facet 

cpright (i) =0 . 0D0 

cplef t (i ) =0 . 0D0 

cptop ( i) = 0 . 0D0 
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Initializing time, opening data files, 



num=12 ! number of system states 

time=o.ODO ! initial time (problem start time), sec 

call sub (xdot, x, time) ! computes dependent variables for initial output 
tprint=time ! initializing time to print output data set 

ptol=.001*delt ! tolerance on time for print 
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VARIABLE. THE ALGORITHM IS DESCRIBED BY THE FOLLOWING 
VECTOR RELATIONS: 


F2=FD (X+H/2, Y+Fl* H/2) , YT=YT+F2*H/3 
F3 = FD (X + H/2, Y + F2* H/2 ) , YT==YT+F3* H/ 3 
F4=FD(X+H,Y+F3*H) , Y=YT+F4*H/6 
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X=X0+H 

CALL SUB (FDERIV, Y, X) 

DO 600 1=1, N 

600 Y ( I ) =YTEMP ( I ) +FDERIV ( I ) * H6 





X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 

Q 


X 

X 

X 

X 

X 

X 

W 


X 

X 

X 

X 

X 

X 

Eh 


X 

X 

X 

X 

X 

X 

DQ 


X 

X 

X 

X 

X 

X 

2 


X 

X 

X 

X 

X 

X 

Cm 


X 

X 

X 

X 

X 

X 

S 


X 

X 

X 

X 

X 

X 

o 


X 

X 

X 

X 

X 

X 

o 


X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 

2; 


X 

X 

X 

X 

X 

X 

w 


X 

X 

X 

X 

X 

X 

w 


X 

X 

X 

X 

X 

X 

PQ 


X 

X 

X 

X 

X 

X 



X 

X 

X 

X 

X 

X 

00 


X 

X 

X 

X 

X 

X 

c 

Q 

X 

X 

X 

X 

X 

X 

PC 

2 

X 

X 

X 

X 

X 

X 


< 

X 

X 

X 

X 

X 

X 

2 


X 

X 

X 

X 

X 

X 

tC 


X 

X 

X 

X 

X 

X 

Eh 



X 

X 

X 

X 

X 

X 

CO 

PC 

X 

X 

X 

X 

X 

X 


+ 

X 

X 

X 

X 

X 

X 

2; 

X 

X 

X 

X 

X 

X 

X 

o 


X 

X 

X 

X 

X 

X 

M 

X • 

X 

X 

X 

X 

X 

X 



Eh 


zc 


X 

X 

X 

X 

X 

X 


< 

>M 

+ 


X 

X 

X 

X 

X 

X 


pz 

PQ 

X 


X 

X 

X 

X 

X 

X 


o 



X 

X 

X 

X 

X 

X 

* 

u 

a 

X 


X 

X 

X 

X 

X 

X 

* 

Eh 

w 

PQ 


X 

X 

X 

X 

X 

X 

* 

2; 

u 



X 

X 

X 

X 

X 

X 


1 — 1 

< 

Q 


X 

X 

X 

X 

X 

X 

00 


X 

w 


X 

X 

X 

X 

X 

X 

00 

w 

Pm 

u 


X 

X 

X 

X 

X 

X 

< 

2; 

u 

< 


X 

X 

X 

X 

X 

X 

CM 

o 

oc 

2! 


X 

X 

X 

X 

X 

X 



Pm 


X 

X 

X 

X 

X 

X 

zc 

Eh 

2; 

PQ 


X 

X 

X 

X 

X 

X 

Eh 

2 

w 

oc 


X 

X 

X 

X 

X 

X 

0C 

H 

u 



X 

X 

X 

X 

X 

X 

D 

o 

PQ 

2 ; 


X 

X 

X 

X 

X 

X 

o 

D-i 


w 


X 

X 

X 

X 

X 

X 

2 


CO 

PQ 


X 

X 

X 

X 

X 

X 


00 

< 

PQ 


X 

X 

X 

X 

X 

X 

Ua 

M 

rn 


2; 

X 

X 

X 

X 

X 

X 

o 

zc 


00 

X 

X 

X 

X 

X 

X 

X 

Eh 


< 

D 

X 

X 

X 

X 

X 

X 

Q 


X 

ZC 

E- Q 

X 

X 

X 

X 

X 

X 

2 

Eh 



PQ 2 

X 

X 

X 

X 

X 

X 

w 

<c 

X 

X 

oc w 

X 

X 

X 

X 

X 

X 





X 

X 

X 

X 

X 

X 

* 





X 

X 

X 

X 

X 

X 

* 





X 

X 

X 

X 

X 

X 

-X 





X 

X 

X 

X 

X 

X 






X 

X 

X 

X 

X 

X 

u 

c 

c 

o 

a 

u 

X 

X 

X 

X 

X 

X 


4-1 

CO 

>1 

to 

X 

P 

c 

•H 

X 

2 

CO 

o 

p 


I — I 

rd 


o 



p 


c 


<15 


TO 


c 


2 


TO 


0) 


p 


>1 


r— 1 


0) 

00 

-p 


03 

£ 

£ 

p 


0 

p 

c 


1 

2 

1 — 1 
CL) 

x 

> 


-p 

„ ^ 

4-1 

no 

"H 


4-1 


0) T5 


__ K 

- 



00 










(D 

p 

4-1 


00 

■ — 










TO 

c 




c 










TO 

-H 

2 


B 

! 










< 

p 

.X 


p 

o 











a 

p 


o 

p 










— • 

p 

r 


c 

0) 











*. 




1 

05 











p 

B 

N 

TO 

l 




00 







c 

05 

1 

03 

p 











*1—1 

p 

o 

p 





MD 






X 

p 

C r> 

V 






' — ' 






1 

a 

o 

X 


00 




p 






o 

r — | 

p 

1 

00 

— 




o 






p 

<d 

Oa 

05 

' — 

X 




p 






a> 

TO 



TO 

l 




CD 



TO 


X 

03 


c 


03 

o 



oo 

u 



TO 


p 

1 

X 

•H 

c 

p 

p 




05 

oo 


03 


p 

p 


03 

o 


0) 



p 

p 





03 

*• 

a) 

B 

"fH 


03 



c 

s. 

00 


x~~i 


a> 

o 

B 


CO 

CM 

1 



o 

V 

— ' 


03 


p 

p 

*H 

p 

• H 

1 — t 

p 



o 

00 

p 


p 



a) 

p 

0 

o 

— 


00 


1 


03 

Q) 

0) 


to 

05 


4-1 

a> 

> 

, 

• — • 


p 

MO 

£ 

£ 

X 

p 

03 

1 

p 


p 

o 

00 

XI 

00 


- — ' 

03 

tn 

p 

TO 

£ 

p 

2 

03 

Pm 

X 

— - 

1 


„ — . 

p 


\ 


> 

p 

N 

0 

p 



o 

o 

> 

00 

o 


to 

a) 

o 


0) 

1 

03 

CD 

, „ 

p 

p 

03 

— 

p 

00 

p 

> 

\ 

CD 

u 

2 

TO 

1 — 1 

CM 

a) 

a) 

p 

to 

a) 

— 

c 

03 

p 

TO 

p 

X 


X 

1 — 1 

03 

05 

cn 

to 

u 

•H 

03 

CO 

p 

2 

o 

p 

a) 

2 

' — 

1 

1 

1 

1 

03 

X 

p 

03 

Q) 

p 

p 


p 

o 

X 

p 

p 

p 

p 

p 

a 

CO 

p 

> 

*H 

0 

<D 

05 

to 








c 

a) 

c 

p 

p 

c 

p 


c 

c 

c 

c 

c 

c 

a 

0 

X 

o 

> — 1 

Q) 

■H 

CO 

p 

0 

o 

o 

o 

0 

0 

o 

u 

p 

u 

m 

05 

p 


•rH 

*H 

-H 

-H 


■H. 

*p 

-H 

\ 

\ 

\ 

N 


2 

to 

u 

CO 

CO 

CO 

CO 

CO 

CO 

CO 

c 

c 

c 

a 

c 

o 

a) 

■H 

c 

2 

a 

c 

c 

c 

c 

o 

o 

o 

o 

o 

p 

p 

rH 

CD 

cd 

a> 

(1) 

a) 

(1) 

a) 

£ 

£ 

£ 

£ 

£ 

X 

♦H 

CLa 

B 

B 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

£ 

p 

p 

B 

•p 

•H 

-H 

*H 

"H 

•H 

■|H 

o 

o 

o 

o 

o 

CO 

5 

m H 

TO 

TO 

TO 

TO 

TO 

TO 

TO 

u 

u 

o 

u 

u 


* 

* 

+ 

u u o o 


o 


o 


179 


common /aerotorq/ t_aero_b ! so that this data is current at print time. (6/27/00) 
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xcv (i+6) =phi (i) /cvdr ! limited Euler angles (deg) 
xcv (i+9) =x (i+9) /cvdr ! body rates (deg/sec) 
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Construction began 4/28/2000, by Doug Nixon, TD55. The purpose 
supply the dynamics xdot=F (x, time) for integration by rksys, a 
fourth-order Runge-Kutta routine. The problem is to set up a 
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The system has six degrees of freedom, or twelve states. The System states are 
defined as follows: 
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r(l)=x(l) ! Changing from mathematical to physical notation, 

r (2) =x (2) 
r (3) =x (3) 



v ( 1 ) =x (4 ) 
v(2) =x (5) 
v (3) =x (6) 
phi (1) =x (7) 
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Aerodynamic drag/torque calculations 
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